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Abstract 

Background: The modelling framework is proposed to study protection properties of 
antibodies to neutralize the effects of the plant toxin (ricin). The present study extends 
our previous work by including (i) the model of intracellular transport of toxin to the 
Endoplasmic Reticulum and (ii) the model of the internalised antibodies (when 
antibody is delivered directly into the cytosol). 

Method: Simulation of the receptor-toxin-antibody interaction is implemented by 
solving the systems of PDEs (advection-diffusion models) or ODEs (rate models) for the 
underlying transport coupled with mass-action kinetics. 

Results: As the main application of the enhanced framework we present a 
comparative study of two kinds (external and internalised) of antibodies. This 
comparison is based on calculation of the non-dimensional protection factor using the 
same set of parameters (geometry, binding constants, initial concentrations of species, 
and total initial amount of the antibody). 

Conclusion: This research will provide a framework for consistent evaluation and 
comparison of different types of antibodies for toxicological applications. 

Keywords: Toxin, Antibody, Cell receptor, Intracellular transport 



Background 

The plant toxin ricin made from the seeds of the castor oil plant is highly toxic to mam- 
malian cells. It is one of the deadliest toxins known and is classified as a potential bioterror 
agent for which no treatment is available. Some promising results have been shown 
recently in the immunotherapeutic approach, i.e. application of antibodies to neutralise 
the effects of ricin [1-4]. With the recent progress in bio-engineering, antibodies with high 
affinity have been generated. The development and production of new antibodies is still 
an expensive process that usually includes extensive experimental studies with continu- 
ous experimental refinement of antibodies properties. Evidently, that such a retrospective 
evaluation of different antibodies aiming at selection of the best candidate may become 
very time and resource consuming. 

In order to reduce this experimental burden a simple (but scientifically consistent) mod- 
elling framework has been recently proposed [5-9]. This framework enables extensive 
theoretical optimization studies to increase the protective potential of antibodies before 
proceeding with targeted experimental studies. 
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The mechanism of ricin intracellular transport involves a number of steps each with 
its complex phenomenology which are well-documented (see [1-4] and Refs. therein). 
Ricin consists of an A (RTA) and B (RTB) chain linked by a disulphide bond. RTB binds 
to a cell surface receptor triggering uptake and retrograde transport to the Endoplasmic 
Reticulum (ER). In the ER, the RTA and RTB chains are separated and the RTA is translo- 
cated across the ER membrane into the cytosol. Subsequently RTA reaches ribosome and 
damages the protein production machinery of the cell resulting in the cell death. In this 
context, the toxin concentration in the cytoplasm near ER becomes the critical quantity to 
estimate the toxicological impact of ricin on the cell and evaluate the protective potential 
of the antibody. This was a motivation to introduce a consistent quantitative characteristic 
for antibody comparison (see below). 

For the sake of parametrisation simplicity the coarse-grained modeling framework pro- 
posed in [5-9] seemingly ignores these fine details of toxin binding and internalisation. 
In fact, it is aimed at capturing the complexity of these processes by means of a small 
number of aggregated' rate constants that can be (or have been) evaluated experimen- 
tally or numerically. Such kind of models becomes a conventional tool in pharmacological 
modeling (for example, see [10] and Refs. therein). From the chemical point of view the 
framework is similar to one well-established in electrochemistry where it is used for esti- 
mation of uptake rates of the heavy metal ions from the environment, see [11-13]. A 
practical application of the proposed models involves a numerical (or sometimes ana- 
lytical) solution of a nonlinear system of PDEs (diffusion kinetics) for a given set (or 
range) of antibody parameters (i.e. binding rates, concentration) to infer the effect of these 
parameters on the protective potential of the antibody. 

In the present paper we extend our previous work [5-9] by refining models for intra- 
cellular transport and chemical interaction of species. Motivated by experimental studies 
available in the literature [1-4] and possible toxicological applications we consider two 
scenarios of antibody delivery. In the first scenario (below refer to as Scenario I) the anti- 
body is placed inside the cell between the ER and cell membrane. The toxin initially is 
delivered outside the cell (in the extracellular domain). Then it moves toward the cell 
and interacts with receptors on the cell membrane. Some of toxin penetrates into the 
cell, where it further moves toward the ER (see Figure la) and eventually interacts with 
antibody. 

In the second scenario {Scenario II) the antibody is delivered outside the cell where it 
initially interacts with the toxin partially neutralising it. The free toxin moves toward the 
cell membrane where it interacts with membrane receptors. Some of toxin penetrates into 
the cell where it is subsequently transported toward ER (see Figure lb). 

To model biologically complex phenomenology of species transport in cell we employ 
the following physical mechanisms: (i) diffusion (for species in extracellular space), (ii) 
diffusion and advection (direction transport via microtubule network) for toxin and 
only diffusion for antibody and toxin-antibody complex inside the cell, (iii) mass-action 
kinetics (to form complexes and describe binding processes). 

Although all parts of the model have been implemented with the standard numerical 
algorithm (an implicit finite-difference scheme [14]), an interpretation of its output and 
parameters tuning requires some domain knowledge (solution of nonlinear systems of 
PDEs). In order to reduce dimension of the parameter space of the models and ameliorate 
their operational application we also considered so-called simplified (or reduced) models 
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Figure 1 The schematic diagram of receptor-toxin-antibody system: (a) - Scenario I and 

(b) - Scenario II. A - toxin, O - antibody, conglutinated A and O in circle - toxin-antibody complex, external 

sphere - cell membrane, internal sphere - ER envelope. 



for each scenario. The simplified (sometimes called compartment) models are based on 
the well-mixed' assumption which states that all species have uniform concentration 
across the modeling domain. This assumption eliminates the a necessity to calculate the 
gradient driven fluxes in the models (i.e. diffusion) leading to a significant simplification 
(translation from PDEs to ODEs). By comparing outputs of complete and reduced mod- 
els we can infer and defensively comment on the trade-off between simplicity and fidelity 
in modeling approach to an evaluation of each scenario. 

As the main application of this framework we present a comparative study of external 
(conventional) and internalised antibodies. More specifically, for the same set of parame- 
ters (geometry, binding constants, initial concentration or initial amount of antibody) we 
calculate the non-dimensional parameter 

S(t) = Mrfcp " )lA6 ^ , (1) 

which is a ratio of reductions of toxin concentration in the cytoplasm due to introduc- 
tion of two kinds of antibodies. Here ujit, p n )\AeQi> ujify p n )\AeQ e ls the concentration 
of internalized toxin near ER for the case of internalized (/) and external antibody (e). The 
time evolution of this aggregated parameter 8 (t) enables a consistent comparison of two 
scenarios of antibody treatment (i.e. either delivered extracellularly or into the cytosol). 
For instance, the condition 8 < 1 indicates that internalised antibody performs better 
than conventional antibody, while at 8 > 1 the conventional antibody outperforms. This 
comparison is the main outcome of our study. 

The paper is organised as follows. In Section 'The models' we introduce models for 
toxin-antibody interaction outside and inside the cell. Results are discussed in Section 
'Numerical results'. Conclusions and summarising remarks are presented in Section 
'Conclusions'. 

Notation 

We use the notation of papers [5-7]: 

T and A - the toxin and antibody; 

C = TA - the toxin-antibody complex (nontoxic); 
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p - the distance to the origin; 

S m = {p : p = p m } - the surface of the membrane of a spherical cell; 

S e = {p : p = p e }, p e > p m - the surface of the external sphere (external surface of 

toe); 

Q e = {p : p e (p m , p e )} - the extracellular domain; 

S n = {p : p = p n }, p n < p m - the surface of the spherical envelope of the domain 
occupied by ER; 

Qi = [p : p e (Pm Pm)} - the intracellular domain; 

ro - the concentration of receptors confined to the cell membrane; 

0 - the fraction of the toxin-bound receptors; 

roO - the concentration of the toxin-bound receptors; 

rn(l — 0) - the concentration of the free receptors; 

ua and uq - the concentrations of the antibody and toxin-antibody complex, 
respectively, in the domains Qi or Q e ; 
ut - the toxin concentration in Q e U Qi; 

- the initial concentrations of the toxin and antibody in Q e and Qi, 
respectively; 

kt, ka, kc - the diffusivity of the toxin, antibody, and toxin-antibody complex, 
respectively; 

k\ y k-i - the forward and reverse constants of the toxin-antibody reaction rate; 

I<2 and k-2 - the forward and reverse binding rate constants of the toxin and receptor 

confined on the membrane; 

k - the toxin internalization rate constant from Q e across the membrane into the cell; 

y - the toxin absorption rate constant describing toxin influx into ER; 

v(p), v > 0 - the toxin advective velocity; 

d t = d/dt, A = P _2 ^(P 2 ^) - the Laplace operator. 

The models 

As mentioned above we study two scenarios of antibody delivery. In Scenario I the 
antibody is delivered inside the cell, i.e. in domain Qi (space between the ER and cell mem- 
brane) while in Scenario II it is delivered outside the cell (space between cell membrane 
and external boundary of the extracellular domain, Q e ). In the Scenario I toxin being ini- 
tially in Q e (extracellular domain) moves toward the cell and interacts with receptors on 
the cell membrane. Subsequently some of toxin penetrates into the intracellular domain, 
Qi, where it moves toward ER and interacts with antibody. Toxin competitively reacts 
with antibody in the extracellular domain {Scenario II) or inside the cell {Scenario /). This 
'blocking' reaction results in a reduced toxin concentration on the ER envelope and is the 
main bio-chemical mechanism of toxin neutralization by antibody. 

In order to extend our results from a single cell model to many-cell systems we impose 
the no-toxin flux boundary condition on the external surface S e of the extracellular 
domain Q e (to mimic periodicity of the infinite system). This also accounts for conserva- 
tion of species in the system and enables consistent simulation of depletion effects (see 
[13] for details). 

It is well-known that binding of toxin to antibody and toxin to receptors are reversible 
reactions, that can be described by the equations of mass-action kinetics (see [5-10] and 
Refs. therein). In the context of study of processes in the Scenario I these equations 
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have been modified to include toxin directional transport and diffusion of all species 
in the intracellular domain (by adding appropriate advection and diffusion terms, e.g. 
see [15,16]). To describe extracellular transport and toxin binding by cell receptors we 
model the cell as a partially absorbing sphere with distributed binding sites that follow 
the Langmuir adsorption model. In the intracellular domain of the Scenario I the toxin 
transport is modelled by the advection-diffusion equation in which an advective velocity 
(drift) is introduced to account for directional transport via microtubule network. For the 
transport of the antibody and antibody- toxin complexes we still use the pure diffusion 
mechanism. We also assume that antibody and toxin- antibody complexes do not inter- 
nalize and never cross the cell membrane or the ER envelope (internal concentric sphere). 
The schematic diagram of receptor- toxin-antibody system is presented in Figure 1. 

Since the process of receptor binding is very rapid (seconds) and toxin internalization 
is relatively slow (hours) we can effectively decoupled the model for extracellular and 
intracellular domains and solve them sequentially. The relative values of time scales for 
receptor binding and toxin internalization have also been validated retrospectively and 
found to be consistent with the initial assumption of domains decoupling (see below). 

The Advection-Diffusion Model used to evaluate Scenario I is an extension of the mod- 
els previously published [4-9]. The following new processes have been included in our 
simulations: (i) the intracellular traffic that captures the diffusion of all species and the 
main features of the toxin microtubule transport, (ii) chemical interaction of all species 
inside the cell (see [15] and Refs. therein), and (iii) process of the toxin transport to ER. 

In order to model toxin transport towards the ER a special boundary condition was 
introduced on the ER envelope. The radius of that spherical envelope was about the size 
of the domain occupied by the ER. It is assumed that the mixed (or radiation) bound- 
ary condition for species concentrations can capture (at least phenomenologically) the 
complexity of toxin kinetics inside the ER [17]. 

Now we can formulate the equations for the Advection-Diffusion Model (refer to as the 
ADM1 model). In extracellular space Q e we have the following set of equations: 

d t U T = K T AU T , P e (p m , Pe), t > 0, 
d p U T = 0, p = p e} t > 0, 

r 0 (2) 
dpUT = — (£2(1 - 0)ut - k-26), p = p m , t > 0, 

Kj 

UT\t=0 =U° t , pG (p m ,Pe), 

d t 0 = k 2 (l - 0)u T - k- 2 0 - kO, p = p m , t > 0, 

(3) 

0\t=o =0, p = p m . 

Function 0 determined from solution by these equations is used as the boundary 
condition in the equations in intracellular domain Qf. 

d t ur = —k\ujUA + k-\uc + ktAut + dpiyur) + 2vujI ' p, 

P e (p n , p m ), t > 0, 

' KT^pUT = rokO — VUj, p = p m > t > 0, (4) 
K T d p U T = (V - V)Uj, p = p n , t > 0, 
UT\t=0 =0, pe (Pn,Pm), 
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d t UA = -k\ujUA + k-\uc + kaAua, P e (p n , p m )> t > 0, 
d p u A = 0, p = p m , t > 0, 

(5) 

d p U A =0, p = p nt t > 0, 
UA\t=0 = P G (Pn>Pm)> 

d t uc = k\u T UA - k-\uc + kc^uc, p e (p n > p m )> t > 0, 

d p Uc = 0, p = p m , t > 0, 

(6) 

a p ^ c = 0, p = p n , t> 0, 
Uc\t=0 =0, pe (Pn,Pm)- 

The set of Eqs. (2)-(6) composes the ADM1 model The main parameter of interest for 
the Scenario I is the antibody protection factor (a relative reduction of the internalized 
toxin due to application of antibody) defined by the expression 

uj(t, p n ;0) 

By definition, 0 < < 1 with the lower values of corresponding to the more 

profound therapeutic effect of antibody treatment. 

To simplify the ADM1 model we employ the well-mixed' assumption [5]. This assump- 
tion implies that all species (toxin, antibody, and toxin-antibody complex) are uniformly 
distributed in the calculation domain for all time, so there is no spatial gradients of con- 
centrations. Under this condition all flux terms disappear from the equations of ADM1. 
The process of toxin internalization (i.e. flux of toxin through the cell surface S m and ER 
envelope S n ) can be modelled as an appropriate rate of toxin removal from Q e and Qi, 
respectively. 

It is worth emphasizing that a reduction to the well-mixed' (or compartment) mod- 
els requires fulfilment of some condition that can be formulated in terms of a smallness 
of the ratio e = t k /tr, where z K and tr is a diffusion and reaction time scale in the 
system, respectively (diffusion can quickly restore any spatial inhomogeneity of species 
distribution maintaining uniform concentration). In the simplest case r K ~ a 2 Ik and 
tr ~ 1/k, where a is a characteristic length scale, k is a scale of diffusivity and k is the scale 
of reaction rate. This is condition for what is called diffusion- or reaction-dominated 
regimes in diffusion kinetics [18]. Indeed, this condition does not always hold and we do 
not assume that it satisfies automatically in the context of our study. Since our models 
are characterized by a number of spatial scales (extracellular and intracellular domains, 
ER envelope), hierarchy of reaction rates and diffusion coefficients, it is very difficult 
(or even impossible) to formulate and validate any general criteria for the feasibility of 
the well-mixed assumption. For this reason we take a heuristic approach in which this 
assumption is validated retrospectively by comparing the output of the models, viz. the 
full (advection-diffusion) and reduced (compartment) models. 

Under the well-mixed assumption the ADM1 model reduces to two ODEs. 

iij = -I<3r 0 (k2(l - 0)ut - k-20), t > 0, 
ur\t=o =4' 

(8) 

0 = k 2 (l - 0)u T - k- 2 0 - kO, t > 0, 



\t-- 



,o=0, 
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(9) 



for extracellular domain Q e and three ODEs for inracelluar domain: 

iij = —k\ujUA + k-\uc + kk^r$6 — ksyuj, t > 0, 
ur\t=o = 0, 

ua = —k\ujUA + k-\uc, t > 0, 
uaU=o = U A' 

uq = k\ujUA — k-iuc, t > 0, 
uc\t=o = 0, 

With/C 3 = 3p2 /( p 3 _ p 3 }j /C4 = 3p 2 /(p 3 _ p 3 }> fe = 3p 2 /(p 3 _ p 3 }< 

We call systems (8), (9) as the Well-Mixed Model 1 (WM1). 

Analogously we formulate two models used to evaluate antibody protection properties 
in Scenario II. For the extracellular domain Q e we employ the Diffusion Model proposed 
in [5-7]: 

d t UT — —k\UjUA + k-\UQ + KtAut, P G (p m > Pe)> t > 0, 
dpUr = 0, p = p e , t > 0, 

KrdpUr = rn{/<2(l - 0)ut k^O}, p = p m , t > 0, 
utU=o = u° T , p e (p m ,p e ), 



(10) 



3 t UA = -k\ujUA + k-\uc + ka^ua, p e (p m , p e ), t > 0, 

d p U A =0, p = p e) t > 0, 

a p ^ = o, p = p m , £ > o, 

WA|t=0 = «a» ^ G (Pm>Pe), 

d t Uc = k\U T UA ~ k-iUQ + KC^Uc, P £ (p m , p e ), t > 0, 
9 p ^c = 0, p = p e , t > 0, 

a p ^c = o, p = p m , t > o, 

Uc\t=0 =0, pe (pm,Pe), 



(ID 



(12) 



(13) 



9,0 = k 2 (l - 0)u T - k- 2 0 - kO, p = p m , t > 0, 
0\ t =o = 0, p = p m , 

with toxin concentration in £2; being described by the Advection-Diffusion equation 

a^^r = ktAut + dpiyur) + 2vut/p, p e (p w , p m ), £ > 0, 
KjdpUj = rokO — vut, p = p m , t > 0, 

K T d p U T = (V ~ V)U T , p = p n , t > 0, 

ur\t=o =0, p e (Pn>P m )- 

We refer to systems (10), (13) as the Advection-Diffusion Model 2 (ADM2). The main 
parameter of interest evaluated with this model is the antibody protection factor 

= u T (t,p n ;u° A )\ Ae n^ 

UT{t, Pn>0) 
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Similarly to the Scenario I we can apply the well-mixed assumption to simplify this 
model and reduce it to a set of ODEs 

iij = —k\ujUA + k-\uc — /<3rn(/<2(l — 0)ut — k-20), t > 0, 
ur\t=o = u° T , 

ua = —k\ujUA + k-iuc, t > 0, 

U A \t=0 = U A , 

(16) 

uq — k\ujUA — k-\uc, t > 0, 
uc\t=o = 0, 

0 = k 2 (l - 0)u T - k- 2 0 - kO, t > 0, 
0\t=o=0 

inQf> and 

uj = kk^roO — yksUT, t > 0, 

(17) 

ur\t=o = 0 

in Qi. We call system (16) and (17) as the Well-Mixed Model 2 (WM2). 

To compare the antibody protection capability in both scenarios (i.e. to infer which 
scenario of antibody delivery is more beneficial for a particular type of antibody) we use 
function 8(t) defined in Eq. (1) in which ujit, p n )\AeQi and ujit, p n )\AeQ e are the toxin 
concentrations near ER determined for the Scenario I and 77, respectively. It is evident 
that the following identity holds, 8(t) = \i\ (t)//ji2(t). 

It is convenient to translate Eqs. (1)-(17) to the non-dimensional form by using scales 
of r* (time), / (length), and u* (concentration). By substituting new variables, p = Ip, 
t = ro = lu*ro, Uj = u*ut> u a — u*ua, uq = u*uc> — u*u^r, i/^ = u*u^, 
k\ = r*uj<i, k 2 = r*u*k 2 , k-i = r*k-i, k-2 = r*k-2, k = r*k, kj = t*ktI~ 2 , y = h* l y, 
v = h~ l v, ka = t*kaI~ 2 , kc = t*kcI~ 2 , h = lh, h = Ik^, h = lh> into (1)-(17) we can 
deduce the same systems, but now in non-dimensional form. Therefore, for simplicity in 
what follows, we treat system (1)-(17) as non-dimensional. 



Numerical results 

We treated Eqs. (2)-(6), (8), (9), (10)-(14), and (16), (17) numerically for t > 0 in spher- 
ically symmetric domains p e (p m , p e ) and p e (p n , p m ), respectively. We solve PDEs, 
Eqs. (2)-(6) and Eqs. (10)-(14), by an implicit finite-difference scheme [14]. To solve 
ODEs, Eqs. (8), (9) and Eqs. (16), (17), we apply the Runge-Kutta scheme. 

Our selection of the values of parameters for the models (2)-(6), (8), (9), (10)-(14), 
and (16), (17) was motivated by the values available in the literature [2,5,19-21] with an 
extended range to allow exploration and illustration of the various transport regimes that 
are possible inside the cell. The following values were used in simulations [18]: u* = 

6.02 • 10 13 cm" 3 , r* = 1 s, r 0 = 1.6 • 10 4 /Sm> where 1.6 • 10 4 is the total number of 
receptors of the cell, / = 10 -2 cm, S m — 47rp^ = 4tt • 10 -6 cm 2 , rn = 2.115 • 10 -3 . Values 
of the dimensionless parameters were the following: p n = 0.02, k = 3.3 • 10 -5 , k\ = 

1.3 • 10~ 2 , k-i = 1.4 • 10~ 4 , k 2 = 1.25 • 10~ 2 , k-2 = 5.2 • 10~ 2 , k a = 10~ 2 , k c = 
10 -2 , iPj — 0, ro = 0.21 • 10 -2 , p m = 0.1, p e = 0.18. Values of kt and y are given in the 
legends of plots. The standard value of v is 0.001. Otherwise it is given in legends of plots. 
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We expect that the chosen values of the parameters were representative enough to illus- 
trate a rich variety of possible scenarios of the evolution of the Receptor-Toxin- Antibody 
system and provide a reasonable estimate of time scales of the associated dynamics. 

The Advection-Diffusion Models enables the fine spatial-temporal resolution of species 
concentration and provide valuable insights into the phenomenology of species interac- 
tion in our system. We found them to be useful tools for understanding the evolution 
of aggregated parameters /xi, /X2, 8 predicted by our simulations. More specifically, by 
means of the Advection-Diffusion Models we were able to identify a number of important 
regimes of toxin penetration that may have significant implications for the assessment for 
the antibody protection properties. 

These regimes can be best described by using the 'transport phenomena terminol- 
ogy which is well-established in chemical engineering, for instance see [22,23]. Initially 
all species are concentrated in spherical layers between the external sphere, cell mem- 
brane, and ER envelope (i.e. (p m , p e ) and (p n , p m )). This layer structure results in strong 
radial gradient of toxin concentration which facilitates the development of diffusion fluxes 
across the system (since cell membrane and ER envelope are penetrable'). Without anti- 
body toxin would eventually be flushed out of the system when it reaches the ER envelope 
driven by the pure diffusion mechanism. The advection velocity will change toxin trans- 
port at one part of the system (from diffusion dominated to advection dominated). How 
much toxin can be blocked by antibody during this transition to ER depends on a number 
of factors. They include the initial reactants concentration (toxin and antibody), relative 
value of toxin advection flux, and reaction time of toxin-antibody binding. It also depends 
on the availability of antibody in the areas of high toxin concentration to maximize the 
effective toxin-antibody binding. The last condition (local availability of antibody for 
binding toxin) is determined by the antibody diffusivity and can be validated by inspection 
of the solutions of Advection-Diffusion Models (spatial-temporal outputs). This valida- 
tion enables justifiable comments on the relative reduction of toxin in the output of the 
system (on ER envelope) and the specific values of antibody protection parameters 
/X2. This approach was employed in our study to identify and investigate the most Vulner- 
able' regimes of toxin penetration and comment on the shapes of functions 112(f). 
Some of these regimes are presented in examples below. 

Furthermore such a chemical engendering' approach provides a clear path for the sys- 
tem optimization. Indeed, from the formal point of view our system (cell) can be cast as a 
chemical reactor with the only output (product), being the toxin concentration at the ER 
envelope. In this context, the aim of antibody treatment is to reduce this output by dis- 
tributing the minimal amount of a given antibody across the system. We anticipate that 
in such settings the optimization problem can be tacked in reduced parameter space by 
using the established framework for chemical reactor design [23,24]. 

The results of the numerical solutions are presented in Figures 2, 3, 4, 5 and 6. The main 
purpose of our study was to estimate the effects of the toxin diffusivity, advective veloc- 
ity, and the 'ER interface' parameter y (rate constant of ER absorption) on the protective 
properties of the antibody in both scenarios. As such, most plots illustrate the effect of 
these parameters. In order to provide a consistent comparison of the simulation outputs 
we apply the same initial concentration of antibody or the same initial amount (but differ- 
ent concentration) of antibody. The equality of the initial amount leads to the following 
relation for the initial concentration 
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„0i „ t ,0e 
U A = K]U A , 



T] 



(pe/pm) 



1 



(18) 



1 - (Pn/Pm) 3 

which follows from the simple geometrical arguments. All calculations were performed 
for p e = 0.18, p m = 0.1, p n = 0.02 for which r\ = 4.871. This condition becomes espe- 
cially important for the comparison of full and reduced models. Unless stated otherwise 
the value of all parameters (except depicted in plots) are assumed to be fixed and equal to 
those given above. 

As mentioned, our simulations reveal a number of interesting regimes of toxin prop- 
agation. In particular, in Scenario I we observed that there was a range of parameters 
(y, kt, v) for which the diffusive flux of toxin, KTd p UT\p=p m = (rokO — vur)\p= Pm > 
becomes negative near the cell membrane, i.e., it is opposite to the advective flux directed 
toward ER. This may reduce the total toxin flux inside the cell and provide more favorable 
condition for toxin-antibody binding. 



(a) 1,0 
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(b) i,o 
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Figure 2 Influence of toxin diffusivity and absorption constant on antibody protection factor. 

(a) Scenario I. Solid and dashed lines correspond to value of uj determined for y = 0.1 and y = 0.05. Toxin 
diffusivity k t : 1 0 -2 (1 ), 1 0 -3 (2), 5 • 1 0 -4 (3). (b) Scenario II. Absorption constant y =0.1, toxin diffusivity k t : 
10- 2 (1), 5 • 10- 4 (2). u° A e = u¥ = 1 in cases (a) and (b). 
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In other words this phenomenon is a direct result of toxin transport being dominated by 
directional advection (due to velocity v). In our simulations this regime occurred for slow 
toxin transport characterized by parameters y, kt> v and the slow forward reaction rate 
for toxin-antibody binding. It is worth mentioning that the similar phenomenon has been 
reported in toxicological studies [17]. We found that the development of this advection 
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Figure 4 Long-time asymptotic behavior of metric parameter 6(f). Solid and dashed lines correspond to 
8 determined for y = 0.1 and y = 0.01 in case u 0 / = u 0 / = 1 , k t = 0.01 . 
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Figure 5 Effect of variation of toxin diffusivity and absorption constant on parameter 8. Toxin 



diffusivity: (a) k t = 1 (T 4 , (b) k t = 10~ 3 . Solid (u 0 / = u 0 / = 1 ) and dashed (u 0 / = 1 , = 4.87) lines 

,0,e _,,0,i _ iN,^^u P / fl 0,e. 



correspond to <5 determined by advection-diffusion models, and bullets (u 0 / = u^ 1 = 1 ) and circles {u^ e = 1 , 



uj = 4.87) to 8 determined by well-mixed models for y :0.1 (1), 0.01 (2). 



dominated' region of toxin transport with reverse diffusive flux has a positive effect on 
antibody protection properties since it effectively increases the efficiency toxin-antibody 
binding (by increasing time available for toxin-antibody reaction). We also solved problem 
(2)-(9) using position-dependent advective velocity v(p) which is zero at the membrane 
and then rapidly tends to a constant value as p -> p n . In this case because of the boundary 
condition (4)2, the diffusive toxin flux near the cell membrane is directed to ER. 

Figure 2a illustrates solution for antibody protection factor for different values of toxin 
diffusivity kj and the ER absorption rate y (with all other parameters being fixed). These 
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results demonstrate the general trend of 8(t) increase with the increase of kt and y. For 
function fi\(t) (Scenario I, Eq. (7)) we found noticeable non-monotonic behavior for all 
kt and y. This implies that there is a time interval for which the protection potential of 
antibody in Scenario I reaches its maximum and then decays. According to Figure 2a this 
time interval is of order of 500 s. 

Figure 2b demonstrates the evolution of antibody protection factor 112(f) (Scenario II, 
Eq. (15)) for a range of parameters kj and y. Although the shape of function fi2 is very 
different (it is always monotonically decreasing) we can conclude that the value of 
drops twice during the first time interval (about 500 s). Simulations for y = 0.1, 0.05 
and kt = 10 -2 , 10 -3 , 10 -4 provide almost identical outputs. This leads to an important 
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conclusion of insignificance of parameters kj and y for the prediction of antibody pro- 
tection factors in Scenario II which is a direct consequence of equation decoupling for 
extracellular and intracellular domain discussed above. 

The plots in Figure 3 illustrate the effect of variation of parameters kj and y on the 
comparative protective skills of the same antibody. This effect is described in terms of 
metric factor 8(t), Eq. (1). Our simulations generally demonstrate the decrease of param- 
eter 8(t) as y or kj decrease. The simulations have been performed for the same initial 
amount of antibody (Eq. (18)) as well as for the same initial concentration (results not pre- 
sented here). We observe that depending on value of kt and y there exists a time interval 
in which 8(t) < 1, i.e. the application of the antibody in Scenario I is more efficient than 
in Scenario II. We found that this interval grows as kt or y decreases. Overall we found 
that the condition 8 (t) < 1 is rather sensitive to the ER absorption rate y. This conclusion 
is intuitively clear since by reducing value of y we also reduce the toxin flux through the 
system facilitating more effective antibody- toxin binding which, in turn, leads to more 
favorable output for Scenario L 

Figure 4 demonstrates some universal behaviors of function 8 (t) that we identified for a 
broad range parameters of our study. Our simulations reveal that for some intermediate 
time (t e (2 • 10 4 ; 10 6 ) seconds) function 8{t) may approximately saturates to a constant 
plateau (depending on values of y and k), before it rapidly decays to zero. This is illus- 
trated by two plots presented in Figure 4 (note different time scales). In particular we 
found that for u% = u^, y = 10 -1 , kt = 10 -2 function 8{t) varies very insignificantly 
(less that 10%, between 4.9 and 5.2) for the time interval (t e [2 • 10 4 ; 5 • 10 5 ]) seconds, then 
it decreases and vanishes as t oo. Note, the limit 8(t) 0 occurs at the unpractically 
long times, so it has more methodologically interest. We do not have a clear explanation 
of the specific numerical value of this saturation limit (8(t) ^ 5); this will require further 
investigation. 

Figure 5 shows the comparative results of the performance of the full Advection- 
Diffusion models (ADM1, ADM2) and the reduced models (WM1, WM2). We compare 
estimations of the metric parameter 8{t) provided by these models for different toxin 
diffusivity kt, ER absorbtion rate y, and initial antibody concentration. We found that 
for the high toxin diffusivity (kt > 10 -2 ) the outputs of full model and reduced mod- 
els are almost identical. The differences between the estimations of 8(t) increases as kt 
decreases and for kt < 10 -4 the reduced models fails to provide any meaningful predic- 
tions (Figure 5b). We found that even for kt > 10 -4 the reduced models usually poorly 
handle the non-monotonic behavior of function 8(t) at short times. At the long times we 
observed that for the broad range of parameters the full and reduced models are reason- 
ably aligned and can yield reasonable estimation of 8 as t —> oo (within 20% accuracy), 
Figure 5a. 

The plots in Figure 6 illustrate the effect of the toxin advective velocity, v, on the 
prediction of the protection factor 8(t), Eq. (1). Two different values of toxin advec- 
tive velocity are used for simulations, v = 10 -3 and v = 5 • 10 -4 with y = 0.1. The 
toxin diffusivity was kt = 10 -3 (Figure 6a) and 10 -4 (Figure 6b) and the same val- 
ues of the initial antibody concentration in both scenarios were used. Only two plots in 
these figures correspond to the well-mixed models, since they do not involve advective 
velocity v. In general we observed the growth of 8(t) as advective velocity increases, so 
Scenario II becomes more favorable. This is intuitively clear because with an increase of 
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v the toxin neutralization time decreases and therefore Scenario I becomes less favor- 
able. We also noticed the improved accuracy of reduced models as v increased. So, for 
given kt and /, the difference between 8(t) determined by Advection-Diffusion mod- 
els and Well-Mixed ones is small only if advective velocity, v, is small. But, for small y, 
it dramatically increases as v grows, while for large y it increases but not so dramati- 
cally. This is due to the switching to the advection dominated transport of toxin in the 
system. 

Conclusions 

We present a new modeling framework to evaluate the protective potential of various 
antibodies. This framework is based on the equation for the diffusion- and advection- 
diffusion transport of all species in the extracellular and intracellular domains and 
include mass-action kinetics for toxin, receptor and antibody. The advection term and 
radiation boundary condition on the Endoplasmic Reticulum envelope were used to 
model the toxin transport via the microtubule network. We estimated the protection 
factor of antibody (relative reduction of toxin concentration near the ER) for two sce- 
narios of operational relevance, viz., when the antibody is delivered externally to the 
cell and when it is delivered directly into the cytosol. To provide a consistent com- 
parison of antibody performance we estimated the evolution of the metric factor 8(t), 
which is defined as the ratio of the protections factors of an antibody for each sce- 
narios. Based on this definition the condition 8(f) < 1 over a time span implies 
that the internally introduced antibody has better protective capability during this time 
while in the opposite case (8(t) > 1) the internalized antibody is more effective. In 
this study simulated both cases. Our models reveal that depending on values of some 
parameters (primary toxin diffusivity and ER absorbtion rate) the plot of function 8(f) 
can significantly change its shape while keeping the same short- and long-time limits, 
by undergoing a saturation stage at the intermediate times with a universal value of 
8^5. 

To reduce the computational burden of the Advection-Diffusion Model (a system of 
nonlinear PDEs) we investigated its possible simplification. By employing the well-mixed 
assumption we reduced the Advection-Diffusion Model to the compartment or rate 
model (set of ODEs) which can be easily solved with any numerical solver. For an extended 
range of parameters we found that predictions of the compartment model can provide 
about ±20% accuracy in estimation of antibody protection potential (depending on values 
of some other parameters), but with decreasing toxin diffusivity limit the compartment 
model become a poor predictor and should not be used in this context. This finding can 
provide some appreciation of a possible trade-off between simplification and fidelity in 
toxicological modelling and can be an important criterion for selection of operational 
models. 

We would like to emphasise that the comparative study reported in the present paper 
can be treated only as the first step for the application of the proposed modelling 
framework to the practical pharmacological studies. More extensive simulations (to pro- 
vide statistically viable outcomes), data fitting and an established strategy for parameter 
optimization using toxicological observations will be necessary. 
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